library(plotly)
library(MASS)

set.seed(123456)

#generate multivariate gaussian distribution
gauss <- mvrnorm(100000, mu = c(0,0), 
         Sigma =matrix(c(1, .5, .5, 1), 2))

#estimate values from sample
distrEst <- kde2d(gauss[,1], gauss[,2], n = 500)

#generate the 3D surface plot
plot_ly(data=data.frame(distrEst),x=distrEst$x,
        y=distrEst$y,z=distrEst$z, type = "surface")